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Abstract 

The particle number projection method is formulated for density dependent forces 
and in particular for the finite range Gogny force. Detailed formula for the projected 
energy and its gradient are provided. The problems arising from the neglection of 
any exchange term, which may lead to divergences, are thoroughly discussed and the 
possible inaccuracies estimated. Numerical results for the projection after variation 
method are shown for the nucleus ^^^Er and for the projection before variation 
approach for the nuclei ^^'^''Cr. We also confirm the Coulomb antipairing effect 
found in mean field theories. 
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1 Introduction 



The simplest approach to the many-body system is the mean-field approxima- 
tion (MFA). The most powerful MFA is the Hartree-Fock-Bogoliubov (HFB) 
theory, where particle-hole and particle-particle correlations are taken into ac- 
count at the same foot. The success of this simple approach is based on the 
large variational Hilbert space generated by symmetry breaking wave func- 
tions. For an atomic nucleus and in the most general HFB approach as many 
symmetries as possible, in particular the rotational invariance and the particle 
number symmetries are broken. In these approaches the quantum numbers are 
only conserved on the average, a semiclassical approach to the full quantum 
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requirement to the wave function of being an eigenstate of the symmetry oper- 
ators. It is well known [1] that the plain HFB approach works well in the case of 
a strongly correlated regime and with a large number of nucleons participating 
in the symmetry breaking process. In the case of the rotational invariance, the 
ideal situation is represented by a strongly deformed heavy nucleus where a 
large number of nucleons contribute to the collective phenomenon of deforma- 
tion. In the case of the particle number symmetry, in which we are interested 
in this paper, the situation is far from being satisfactory because only a few 
pairs of nucleons ( four to five) in the vicinity of the Fermi surface are thought 
to cause the phenomenon of nuclear superfluidity. Thus, it is well known that 
the treatment of pairing correlations in the mean field approach produces an 
unrealistic transition from the superfluid to the normal phase along the yrast 
band, not expected in a finite system. In order to get an improved description 
of the nucleus it is clear that one has to go beyond the mean field approx- 
imation to include correlations. This task can be achieved in a simple and 
effective way within the particle number projection (PNP) [2] formalism, an- 
other possibility is to build correlations by the Random-Phase- Approximation 
(RPA) [3]. The latter is, however, from the applications point of view, very 
limited, since the required calculation of the ground state correlations can be 
performed only for very simple models and/or separable interactions. Con- 
cerning this point one has to keep in mind that not only the RPA but also 
other sophisticated theories have been developed only in simple models or 
with separable forces. Thus the first calculation with exact particle number 
projection before the variation was done about twenty years ago [4] with the 
pairing plus quadrupole Hamiltonian and the Baranger-Kumar configuration 
space. 

In the last years in an attempt to achieve some progress in calculations with 
effective forces and large configuration spaces, the Lipkin Nogami method [5] 
has experienced a revival. The great advantage of the Lipkin Nogami method is 
that corrects, in an approximate way, for the zero point energy associated with 
the breakdown of the particle number symmetry at the mean field levels i.e. 
the calculations do not get much involved as compared to the full HFB ones. 
Some applications of this method with effective forces have been done with the 
Skyrme force [6,7], the Gogny force[8] and the relativistic mean field approach 
[9]. The Lipkin Nogami method on the other hand is just a recipe and cannot 
be obtained from a variational principle. It is therefore desirable to develop an 
exact particle number projection which can be applied to effective forces. The 
purpose of this paper is to develop an exact particle number projection before 
(and after) the variation for density dependent forces and in particular for the 
Gogny force. As we shall see, the exchange terms play an important role in 
PNP, in fact the ncglcction of any exchange term may lead to divergences. 
As it is well known, in general, a zero range force is unable to provide good 
pairing properties. As a matter of fact, normally the pairing terms of the force 
are neglected in order to get reasonable Hartree Fock results. To perform 
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HFB calculations a pairing force has been usually added ad hoc. The finite 
range of the density dependent Gogny force, on the other hand, provides good 
pairing properties without neglection of any exchange term. In fact, the force 
was designed with this intention[10] and this property makes the Gogny force 
one of the few effective density dependent forces for which particle number 
projection is feasible without problem due to divergences. All the Skyrme 
parametrizations (with the few exceptions of those which provide reasonable 
pairing correlations) and all relativistic mean field approaches may give in 
general divergent contributions to the projected energy. 

In section 2, we introduce the particle number projection method for non- 
separable density dependent forces. In section 3 we present the formalism for 
the density dependent part of the interaction and discuss some alternative 
prescriptions to this dependence. The gradient of the projected energy to 
be implemented in the calculations is developed in section 4. In section 5, 
we present results for the projection after variation (PAV) method for the 
nucleus ^^^Er and discuss the problems associated with the divergent terms. 
Results for the variation after projection (VAP) method for the nuclei ^^'^^Cr 
are discussed in section 6. 



2 Pcirticle Number Projected Hcirtree— Fock— Bogoliubov Theory. 

The Particle Number Projected Hartree-Fock-Bogoliubov theory is the sim- 
plest symmetry conserving mean field approach which takes into account 
particle-particle correlations. This theory has been derived in the past for 
non-density dependent interactions and a detailed formulation has only been 
given for separable forces like the Pairing Plus Quadrupole model [4] . In order 
to generalize this theory to non-separable and density dependent forces, like 
the Gogny or the Skyrme force, and to introduce the pertinent notation we 
shall also summarize the known formulation. 

Let 1$) be a product wave function of the HFB type, i.e, the vacuum to the 
quasiparticle operators given by the Bogoliubov transformation 

k 

where (ca;,c|,) are the particle annihilation and creation operators (usually in 
the Harmonic Oscillator basis). This transformation does not conserve the 
particle number symmetry and the wave function |$) is not an eigenstate of 
the particle number operator N. One can generate, however, out of 1$), an 
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eigenstate, |^iv), of N by the well known projection technique [1] 



27T 

1*^) = P^|$) = —l d<pe'^^^-''^ 1$), (2) 





where the integration interval can be reduced to [0, vr] whenever the intrin- 
sic wave function has a well defined "number parity" [1] quantum number. 
Furthermore, the integral can be discretized in a sum using the Fomenko ex- 
pansion [11] 

pA^= l^e^.K^-iv)^ ^. = 7/, (3) 
^ 1=1 ^ 



where L is the number of points used in the calculations, for L ^ 00 one gets 
the exact solution, usually L = 8 is used in the calculations. projects only 
on one system of particles, in the general case, however, one has to perform 
a double projection, for the number of neutrons and protons. In this case the 
projected matrix element of any operator O, such that [O, TV] = 0, is given by 
[4] 



where we have introduced 



yir = — , (5) 



and the O-operator overlap 



with 1$) = 1$^)!$,^). Another useful quantity is the norm overlap, which is 
given by 

xi^ = ($,|e^'^'^(^^-^^)|$,), (7) 
clearly, yi^ and xi^ are related by yi^ = y^xi^, with y^ = (Et=i ^ir)'^- 
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As it will be illustrated for the Hamiltonian operator, any operator overlap can 
be calculated using the generalized Wick theorem [12]. This theorem allows to 
express matrix elements of the form exactly in the same way as the 

ordinary Wick theorem, as the sum of all possible contracted products, where 
1$) and 1$) are quasiparticle vacua. Thus, in particular 



- ($|4c,J$)($|4c,3|$) + ($|4clj$)($|c,,c,3|$)] (8) 

Looking at eq. (4) we see that this theorem allows to express projected ex- 
pectation values in terms of mean field matrix elements. In eq. (8) the first 
term is the Hartree term, the second one the Fock term and the last one the 
Bogoliubov term. ^From this expression we see that problems may appear 
when ($1$) ~ 0. This point will be discussed later on and its consequences 
for the effective force will be analyzed in Appendix B. For the PNP method, 
the basic building blocks of Wick's theorem are the (/^-dependent generalized 
density matrix and pairing tensors 



= W4:#> = (^l^SiC^ = (v.(^)v-) (9) 

($1$) ($|e^^^|$) ^ ^'^i^' 

($1$) ($|e*'^^|$) ^ 

where we have used the notation |$ >= e*'^^|$ >, the indices {k, k') belonging 
to the same isospin channel r. The matrices W(</?) and V(</?) are related to U 
and V of eq. (1) by 

U{ip) = U + V*X{ip)* V{cp) = V + U*X{ip)\ (12) 



and the contraction X{(p) is given by 



= ^ (13) 



here again the indices v belong to the same isospin channel and 

T22 {v) = cos{^) -1 - i sin((/?) • N^^ (14) 
T2i(99) = -isin((^) -TV^o (15) 
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with N^^ = WU - V^V and iV^o = jj^y* _ ytf/*. The matrices N^^ and A^^o 
represent the number operator in the quasiparticle basis of eq.(l). 

The nuclear Hamiltonian can be written as 
H^T + Vdi + Vdd 

= E *fclfe2 4iCfe2 + I Yl Vkik2k3kAA2^k4Ck3+yDD (16) 

kik2 kik2k3k4, 

where the density independent part, Vdi, represents the nuclear interaction 
(Coulomb included) with the exception of a possible density-dependent term 
(the well known density dependent term of the Skyrme and Gogny forces) 
which is represented by Vdd- In our case Vdd is given by the last term of 
eq. (A.l). 



2.1 Particle Number Projected Energy for the non-density dependent parts 
of the interaction. 

With the definitions above we can formulate the particle number projected 
theory for non-separable interactions, like the Gogny force (Appendix A), in 
a similar way as in the mean field approximation. Using eq. (4) and the gen- 
eralized Wick theorem, the density independent part of the projected energy 
is given by 

L L . 

+ 1 (r'^'^((^,jp((^,j + r-((^,jp((^,j + r-((^,jp((^,j 

+ r''''i^iM^i.) - Ai°-(^,J«:°i(^,J - Ai°''^(^,J«:°i(^,j)} , (17) 

where the traces run over the configuration space, here assumed to be the 
same for protons and neutrons, and 

^kl'ksi'Plr') = E ^klk2k3k4Pk4k2('Pl,r), (18) 
^k!_k2M = ^ E ^klk2k3kyk3k4M- (19) 
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In T'^'^' the indices (fci, ^3) belong to the isospin r and (/c2, ^4) to r', in A^^'"^ 
all indices belong to the isospin r 0. Since J2t=iyiT = we can write the 
proton or neutron projected density as 



with {k, k') belonging to the isospin r, and the projected energy as 



^ r 1^ = 1 ) 



where 

L 



^fcifcl = E VK'^kik^A'fK') = ^ ^krk2hikiPklk2 ■ (22) 
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and r 7^ r'. By (Pr') in T^'^^' and in some notation to be introduced later 
on ( see also eq. 20 ) we want to indicate that these quantities are ipi^, inde- 
pendent, i.e., the integration on ipi^, has been performed. In these fields, the 
indices {ki,k^) belong to the isospin r. The expressions (20) and (21) can be 
rewritten taking into account that 

yL-i=y!, pi^L-i) = pK^i), ^^''ivL-i) = e-''^^K'\^ir. (23) 
In this way we can simplify the sum in eq. (20) for the projected density to[^ 



[L/2] 

p''^ = Re{y^,Y.Ci.^iyifir)} , (24) 

where we have taken into account that yi_^ = yl ■ xi^. By [L/2] we represent 
the integer part oi L/2 and 



^ Notice that the pairing field A does not mix isospin. 

^ Notice that the sum in eq. (3) from / = 1 to L can be written as a sum from / = 
to L - 1. 
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1 if /t- = L/2 and L even 

2 otherwise. 



(25) 



In the same way, for the projected energy, eq. (21), we obtain 

T lr=0 ^ J 



3 Particle Number Projection for the density dependent part of 
the interaction 



Density dependent interactions have a term which simulates a G-matrix through 
an explicit dependence of the nuclear density. This term, Vdd, was introduced 
in the mean field approximation. In the MFA only diagonal matrix elements 
between product wave functions are needed to calculate the energy and, con- 
sequently, Vdd is constructed to depend on the mean field density. In theories 
beyond mean field, the density dependent part is given by 



4 1 : : (27) 

where [p(r)] indicates the explicit dependence of Vdd on a density p(r) to 
be specified. Looking at these expressions it is not obvious which dependence 
should be used. There are two more or less straightforward prescriptions [13] 
for p(r) : 

(1) The first prescription is inspired by the following consideration : In the 
MFA, the energy is given by (<I>|i/|$)/($|<I>) and Vdd is assumed to de- 
pend on the density ($|/3|$)/($|$). On the other hand, if the wave func- 
tion which describes the nuclear system is the projected wave function 
\^n)-i we have to calculate the matrix element n\Vdd\^ n) / n\^ n) 
(see the middle term in eq. (27)). It seems reasonable, therefore, to use in 
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Vdd the density p{f) — n\p\^ n) / n\^ n) — P^^{r), i-e. the projected 
density. 

(2) The other prescription has been guided by the choice usually done in the 
Generator Coordinate method with density dependent forces [14]. The 
philosophy behind this prescription is the following: to evaluate eq. (27) 
we have to calculate matrix elements between different product wave 
functions |$) and |<l>) (|<l>) = e*'^^|<l>)) (see last term in eq. (27)). Then, 
to calculate matrix elements of the form 

(iM ,28) 

we choose the mixed density 

to be used in Vdd- We shall call this approach the mixed density pre- 
scription. 

Both prescriptions have been tested with the Gogny force in the Lipkin Nogami 
approach [15] and practically no difference was found in the numerical appli- 
cations. One should notice that in the second prescription p(r) depends on 
the angle (p at variance with prescription 1. 

We shall now proceed to evaluate the projected energy with both prescriptions. 

3.1 Prescription 1: Projected Density. 

Using the projected density prescription, the projected energy of the density- 
dependent term {E£]j) is obtained in a simple way. One just has to use in Vdd 
the projected density. In the case of the Gogny force (see Appendix A) which 
we are considering here, the value xq — 1 for the spin exchange parameter 
guarantees that there is no contribution from Vdd neither to the pairing field 
A nor to the Hartree-Fock field of the type F"^"^, i.e., without isospin mixing. 
We find 

= ^r^(f •p^' + f'^'^'-p^'') , (30) 

where 

fS'=E(«..ipn),.«„/rx. (31) 

fe2K4 
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with T ^ t' and 

[vddIp''])^ , , = {kMVDD\p'']\kzk^) - {kMyDD[p'']\kJ^^) , (32) 

the symbol " is used to specify that we are deahng with the DD part of the 
interaction. 



3.2 Prescription 2: Mixed Density. 



We assume in this section that Vdd depends on the mixed density 



(33) 



The contribution of Vdd to the energy is given by eq. (27) with [p(r)] = 
Pip^ip^{r). Since the mixed density depends on the angles, we cannot simplify 
the double integral to a single one as we did before. Replacing the integrals 
by sums and using (23) again, we obtain 



r /[L/2] [L/2] \ 

= ^ E E E yi.A.i^> ■ Tr (f-'[p(vp,j, p(^,^,)] • pi^i^)) 

TT^r' { \lr=0 l,,=0 J 

/[L/2]-l [L/2]-l \ ^ 

+2-Re[ E yir E ylTr(T^^'^P^'^^r),P*i^i^,)]-pi^ir))\} (34) 
where [L/2] represents, as before, the integer part of L/2 and 



1 if = or L/2 and Ir' = or L/2 (L even) 

1 if = and 1^' = (L odd) 

2 otherwise. 



(35) 



In this expression we have furthermore introduced the following notation 



hl'k3[Pi'Plr),P{'Pl,')]^Yl {'"DD [p(<^'J'P(<^'.')])fc^fc^fc3fc/fe4fe2(</'i./) 

^ll'ks[PiVlr),P*iVl,,)]^Yl (vdd [p(<^'J'P*(<^'.')])fc^fe,fc3fc/fe4fe2(<^'.')- (36) 

fe4fc2 
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Here again the indices (^1,^3) belong to the isospin r and (^2,^4) to r' 
(r 7^ r') and (vdd pivir): PiVi 1) I is the density-dependent matrix 

element calculated with the density p'^{(pr,r) + p'^ (ipT-',r) in Vdd, similarly, 
(vdd p{'Pir),P*{'Pi^,) ) is calculated with p^{(pr,r) + {p^' {ipr' ,r))* in Vdd- 

If we compare the expression (34) with the result obtained for E^^ with the 
projected prescription (30), it is easy to observe that (34) implies a more elab- 
orated calculation. With the projected prescription, we only have to calculate 
the projected density ^ the field r[p'^] and then one easily gets E^j-,. But 
in order to evaluate expression (34) it is necessary to calculate [1? jl + 2) 
different fields. For example, for L = 8, we have to calculate 34 fields. 



4 Restoration of the particle number symmetry before the varia- 
tion in HFB theories with non-separable forces. 



Using equation (26), we can perform a restoration of the particle number sym- 
metry, projecting after the variational Hartree-Fock-Bogoliubov equations 
have been solved. In this case the HFB wave function is not self-consistently 
determined. In particular, if the HFB calculation collapses to the pairing un- 
correlated HF one, the PNP method does not provide a better approximation. 
In this section wc shall derive the projected variational equations in order to 
solve the self-consistent problem, that means, the variation after the projec- 
tion method. According to the Ritz variational principle, the intrinsic wave 
function |$) has to be determined by minimization of the projected energy, 
i.e., 

bE^ = b ^ ' . ^ ' = (37) 



The variational Hilbert space generated by product wave functions can be 
parametrized by the Thouless theorem [1], which allows to write any product 
wave function, |$), in the form 

|$)=A^exp(Vc,.at«t)l^o) (38) 



where |$o) is a reference product wave function non-orthogonal to |$) and 
Cp^y are the variational parameters. We shall first study the variation of the 
density independent part of the interaction, eq. (26). An arbitrary variation 
of the energy is given by 
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I E + ^-c-) (39) 



2 



where we have introduced the projected gradient Ti^j^. Here the indices 
(/i, z/) run over the proton and neutron configuration space. The solution of 
the projected variational problem can be obtained by the conjugate gradient 
method [16]. 

Using the generalized Wick theorem [12], H^j^ can be written as 



with the index r corresponding to the isospin of the quasiparticle operators 
(cKju, a^). We have furthermore introduced X^"^ — yiT'^iVir) ^^<i 

+ Wt(<^^jAi°'^((^,jW*((^,J - Vt((^,JA°i'^((^,jV*((^^j)^^ (41) 

with 

M^.J=i + r'''^' + r^(^.J, (42) 

again r ^ t\ and 

+ iTr| (v^ (t + r-'^-yC/ -U^{t + r-^^') \/) (A'(^,J - X''^) |43) 



where 



£^Lr(^J - ^Tr(r-(<^,Jp-(^J - A^°'-(<^J/.°^'-(<^J) (44) 



Taking into account eq. (23), we can simplify the gradient expression to the 
region [0, L/2] in the following way : 
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[L/2] 

nf/r = y^oRe E Cw^ir {K^A^w) + s-{<pw) - x^j) } . (45) 

It=0 

These expressions apply to the non-density dependent part of the interaction. 
To calculate the gradient of the density dependent part of the interaction we 
proceed as with the evaluation of the energy, accordingly we shall distinguish 
the two prescriptions adopted. 



4-1 Prescription 1: Projected Density. 



In a similar way as in the case of the density independent part of the interac- 
tion, the variation of E^j^ is defined by SE^j^ = | ('^^^-'^ ^C**^ + h.c.^. 

Here the expression of is somewhat more complicated due to the explicit 
dependence on the density of the interaction : 



Using eq. (A. 3) one can calculate very easily : 



(47) 



Using the notation SF^f,, for the rearrangement term 

( $ I I $ ) 



where the indices {k, k') and v) in eqs. (47,48) belong to the isospin r, we 
can write the projected gradient of Edd in the following way : 



'n'X = vlR^^ E Cu^ir -Hf^A'Pir) + s-{^iXx,A^ir) - <;) (49) 

lr=0 I ) 

where 
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-Vt(^,j(^^^' + 5p)V(^,j) (50) 



and 



(51) 



This is a similar expression to the projected gradient of the non density- 
dependent part, (45), they differ only by the term E'^{(pi^), which does not 
contribute in this case because the specific value of the coefficient xq in the 
Gogny force makes it zero. 



4.1.1 Prescription 2: Mixed Density. 

We shall now derive the gradient of the density dependent term with the 
second prescription. The projected gradient is given by an expression similar 
to eq. (46) but with Vdd depending on the mixed density, eq. (33), and the 
rearrangement term given by 



As before, we obtain 



[L/2] [L/2] 

[L/2]-l [Z,/2]-l N 

1^=1 l^,=l ) 
- E^udK' (53) 

where 
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'f^fvA'P^r->'-P*i^,) and E^jj{ipi^,ip*^) in eq. (53) are given by the corresponding 
expressions in eq. (54) but with the replacement p{<Pi^,) —>■ {p{(piA)*. 



5 The Projection after variation method and the divergences prob- 
lem. 



In this section we discuss some problems that may arise in PNP calculations 
as well as some theoretical results. As we have seen in eq. (8), the extended 
Wick theorem allows to factorize matrix elements of the two-body part of the 
interaction in terms of the generalized matrix density and pairing tensors. The 
matrix element of eq. (8) takes a simple form in the canonical basis. Using 
eqs. (9-11), we obtain 



(^1 4i «L«fc4«fc3l*) = Pk3krPkik2 - Pk^krPksk^ + I^llk2'^llki 



4. 



+ 



Uk-.Vk-.e 



2 

k2 

2iip 



kzk\"k2ki 



}k4k\"k3 k2 



(55) 



The dangerous cases are for (p — 7r/2 and ul — when the denominators 
become zero. Taking into account that ($|$) = nfc>o('"i + "^^6^"^)! easy to 
sec that as long as there is only one pole or two poles but ki ^ k2 { or ^ ki ) 
there are no problems because they will be cancel out by the multiplying 
overlap ($|$). The only case where problems arise is with the matrix elements 
of the form ($|a|.a^a^afc|$), because in this case no cancellation with the 
norm overlap ($|^) is possible. It is easy to check, however, that in this case, 
though each term of eq.(55) by itself is divergent, their sum gives a finite 
contribution, namely vl ■ e^"^ ■ Yim>o,rny^ki'^m + "^m^^"^)) which behaves properly 
when vl = 1/2 and (p = 7r/2. That means, as long as all three terms of the Wick 
factorization of eq. (8), i.e., the Hartree, the Fock and the Bogohubov term, 
are taken into account no divergences appear in the PNP formalism. Notice 
that the arguments given here are quite general and independent of the kind of 
force (Coulomb force included) used in the calculations -see Appendix B for a 
detailed discussion. This point has already been noticed in ref. [17]. Obviously, 
one body operators will never cause problems to the PNP. 
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In the case of the density dependent forces, one must be aware that this 
dependence itself may cause problems. As we have seen before, two different 
prescriptions have mainly been used to treat the density dependence of the 
Hamiltonian. As discussed in Appendix B, the projected density prescription 
does not have divergences and the mixed density prescription has an integrable 
divergence. For the moment, we use the projected density prescription (PAVl), 
at the end of this section we show results with the mixed density prescription 
(PAV2). 

Unfortunately most of the calculations performed so far with effective forces 
neglect exchange terms. Some of them because different forces are used in 
the particle-hole and the particle-particle channel (like most of the Skyrme 
force parametrizations and the relativistic mean field calculations). The con- 
tributions to the p-p interaction from the force used in the p-h channel being 
neglected and vice- versa. Others, though with the same force in both channels, 
like the Gogny force, neglect some exchange terms (Coulomb among others) 
just to save CPU time. Only recently HFB calculations including all terms in 
a triaxial basis has been performed with the Gogny force, see [18]. 

Obviously, the most interesting issue to address is to investigate how large the 
divergences are in the case that some exchange terms are neglected. Though 
this point will ultimately depend on the particular force used in the calcu- 
lation, we shall illustrate it here for the Gogny force where we are able to 
perform the usual calculations which neglect some exchange terms and the 
exact ones. To be precise, we present two types of calculations, the first one 
is the usual one used in calculations with this force in a triaxial basis [19,20], 
namely, the pairing field is calculated only from the Brink-Boeker term, i.e. 
all additional pairing exchange terms are neglected and the Fock term of the 
Coulomb potential is either neglected or calculated in an approximate way. 
We shall call this approximation HFBs, s for standard calculation -see ref. [18] 
for further details. The wave function solution of the selfconsistent HFB equa- 
tions with this approximate Hamiltonian will be denoted $5. In the second 
calculation, the exact HFB equations (HFBe) with the Gogny force are solved, 
i.e., al/ exchange terms are considered in the variational principle -see [18] for 
further details. The wave function resulting from the selfconsistent solution 
of the HFBe equations will be denoted $£. In a later stage the intrinsic wave 
functions obtained in these ways are projected onto good particle number in 
order to calculate the projected energy. 

In the calculations, we have expanded the quasiparticle operators in a triaxial 
harmonic oscillator basis. The configuration space was chosen by the condition 

a^Ux + ayUy + azUz < No (56) 
where the coefficients Oj depend on the relations between the axes q — R^/ Rx 
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and p = Ry/Rx in the form = (gp) , ay = q^/^p'"^/^ and = p^/^q~'^/^. 
This basis has been symmetrized with respect to the simplex operation and is 
an eigenstate of the operator TI2T with 112 = Pe"*'^'^"', T is the time-reversal 
operator, and P the parity operator [21]. In the calculations we take Nq = 11.1. 
To take into account the rotational motion we add to the Hamiltonian the 
cranking term —cuJ^, as in ref. [20,18]. 

The most stringent test one can perform is to study the yrast band of a nu- 
cleus, ^^'^Er for example, since the pairing correlations and the quasiparticle 
occupancies Vk strongly depend on the angular momentum. To separate the 
different issues we shall first study only the convergence of the projected en- 
ergy in terms of the parameter L of the Fomenko expansion. In order not to 
mix this issue with the divergence problem associated with the neglection of 
the exchange terms, we shall solve first the exact problem, i.e. the HFBe case. 
To illustrate more clearly our results we shall perform particle number projec- 
tion after the variation (PAV), i.e., first the cranked HFBg equations with the 
Gogny force are solved, as in ref. [18]. This provides us with the wave function 
$e- Then a projection on the particle number (PAVe) is performed onto the 
wave function $e {ioi different L values ) to calculate the projected energy. In 
fig. 1, we display the results of the transition gamma ray energy in the PAVe 
approach for different L-values. We find that the points for the L-dependent 
calculation are on top of each other, which indicates a clean convergence and 
that L = 8 is already a good approximation to the exact PNP. This value is in 
agreement with earlier calculations [4,22]. We find, by the way, a good agree- 
ment with the experiment [23] up to / = 18^ and a considerable improvement 
to the HFBg approach (dashed fine). For higher spin values, the agreement 
with the experiment is not as good because the proton pairing energy has 
already collapsed to zero. 

We now turn to the problem of the divergences in calculations where the 
exchange terms are neglected. As we have seen in eq. (55) and in Appendix B, 
the poles appear whenever -u^ = f| = 1/2 and Lpi = tt/2. The first condition, 
obviously, can not be eluded. In the non-selfconsistent PAV method because 
the occupancies vl of the wave function are already determined by the solution 
of the HFB equations and in the self-consistent VAP because it would be 
against the variational principle. To investigate the behavior of the occupancies 
vl, we have solved the sclfconsistent cranked HFBs equations along the Yrast 
band for the nucleus ^^^Er. The wave functions $5 determined in this way have 
been analyzed in the canonical basis. In Fig. 2, we show the smallest values 
of the quantity ul — vl, for protons and neutrons, found in as a function 
of the angular momentum. At small spins and up to spin 8 h, the quantity 
u^ — vl for neutrons is rather small (around —0.05), at spin 10 h it reaches the 
critical value ul = vl and from this point on it steadily increases up to the 
maximal value at a spin value of 28 h. This behavior is typical of the high spin 
regime. At small angular momentum, the system is well paired and we may 
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Fig. 1. E^[MeV] versus angular momentum calculated for different values of L in 
the Fomenko expansion, together with the experimental values and the HFBg ones 
(dashed line). Notice that all points for different L values lie on top of each other 
for all angular momenta. 

find a pair of nucleons close to the Fermi surface with m| v\. As the spin 
increases the pairs align (states start being predominantly occupied or empty) 
and at high spins, when the pairing collapse takes place, the states become 
either fully occupied or empty. For protons the situation is even worse, we find 
two critical points, one around / = VlTi and the other around / = 22^. 



> 
I 




8 12 16 20 
Angular Momentum ( R) 

Fig. 2. The smallest values of (^1 — u^) along the Yrast Hne of ^^^Er for the HFBs 
in the canonical basis. 

The second condition for having the pole, i.e. that = 7r/2, is in principle 
avoidable, since ipi — t:1/ taking L odd in the Fomenko expansion should be 
enough. This recipe, as we shall see below, helps in some cases but it does not 
guarantees that one is error free. The reason is very easy, to find convergence 
in the Fomenko expansion one must take L large enough, but then one may get 
too close to the critical 7r/2. That means in this way one avoids to get nonsense 
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but the results may not be reliable. As we have seen above and in Appendix B, 
in the neighborhood of the poles we will get separately big contributions to 
the Hartree-Fock energy and to the pairing energy. If all exchange terms have 
been taken into account in the calculations, the sum of the divergent terms 
provides a finite value. This cancellation, however, is not guaranteed if some 
exchange terms have been neglected, as in the HFBg approach, which we now 
discuss. In Table 1, we show for the nucleus ^^"^Er the differences in the pairing, 
eq. (B.14), Hartree-Fock, eq. (B.12), and total, eq. (B.l), projected energies 
calculated for L = 8 and L = 20 (column labeled L-EVEN) and for L = 9 
and L = 21 (column labeled L-ODD) as a function of the angular momentum. 
In the even-L calculations in the PAVg approach, wc find large AE values 
in the pairing and the HF energies around I = 10, 12^ and I = 24h. These 
deviations do not cancel each other and give rise to large AEtot, indicating 
that no convergence is reached. For odd-L values the situation is much better. 
The largest AE'^t is around 20 keV which is more than one order of magnitude 
smaller than in the even-L case. Lastly, in the last column the difference of the 
total energy between L = 8 and L = 9 is displayed. We clearly observe that 
again large values around / = 10, 12h and / = 24/i are obtained. Interestingly, 
a comparison with fig. 2 shows that these spin values are those where u\ — v1 is 
very small, i.e. the poles. The neutron system is responsible for the increase at 
/ = lOfi, see second column, and the proton one around / = V2h and / = 2Ah. 

lYiom. this table, we conclude that even-L values should be avoided and that 
the total energy for odd-L is a function of L. Usually, we are not interested in 
total energies but in relative ones, for instance in the 7-ray energy. Now the 
question is whether we can find a kind of plateau condition that provides us the 
optimal L- value for all /-values. In Fig. 3 we display the quantity Ei^ {I) / E}:^ {I) 
as a function of L. Wc have chosen the value L = 11 to normalize the 7- 
ray energies because it is large enough that in the absence of poles a good 
convergence would be found, and small enough that in the presence of poles 
we do not come too close to the critical ipL = ti/2. We found a good plateau 
for spin values 2 to Qh and a rather good one for I = 16, 18 and 20^, the worst 
cases are for / = lAh and / = 24h. This behavior is again strongly correlated 
with the nodes in Fig. 2. We observe that the largest deviations of the most 
reasonable values are about 7 %. It is clear that more sensitive quantities like 
second moments of inertia will show larger uncertainties. 

All calculations of projected energies shown up to this point have been done 
with the prescription of eq. (30) for the density dependent part of the force, 
namely the projected density. 

We shall now discuss the convergence of our results with respect to the density 
dependence of the Hamiltonian used in the calculations. In order to isolate this 
issue, we do not neglect any exchange terms in the calculations, that means we 
solve the HFBg equations. This provides us with the wave function $e which 
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Table 1 

Differences (in keV) between tlie projected energies calculated with L = 9 (E'f ) and 
L = 21 (^21) and L = 8 {Ef) and L = 20 {Ef^) in the Fomcnko expansion for 
the nucleus ^^^Er, calculated in the PAVs approach. The last column displays the 
difference between the total energy calculated with L = 8 and L = 9. 

we use to calculate the projected energy, which can calculated according to 
two prescriptions. Notice that $e is different from $5 and that the levels 
occupancies of both wave functions are not necessarily the same. As it is 
shown in section B.1.1 the prescription using the projected density, which 
we shall label DDI, was divergence free, while the second prescription, the 
mixed density one (section B.1.2), from now on called DD2, was showing an 
integrable divergence. As in Table 1 we shall calculate separately the pairing 
and the Hartree-Fock contributions to the total energy. Since the density- 
dependent term does not contribute to the pairing energy (the parameter xq 
in the Gogny force is equal to 1) this energy is independent of the prescription 
used for the density dependent term. In Table 2 we display the differences 
/\^8-20 q£ projected energies calculated with the pairing energy, eq. (B.14), 
the HP energy, eq. (B.12) and the sum of both, the total energy. At low 
and medium spins, the difference of pairing energies for L = 8 and 20 is not 
equal to zero indicating the known divergences. We find the largest value for 
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Fig. 3. The transition energies E-y[MeV] for a given angular momentum as a function 
of the parameter L of the Fomenko expansion. 

I — 4h indicating that the wave function $e has a level occupancy at this 
spin value such that ~ v^. Inspection of the canonical basis indicates that 
the orbital is a neutron one, the proton system in this case does not have 
orbitals with vl close to 1/2. At high spins, the wave function <l>e is not paired 
in the proton system (see ref. [18]) and in the neutron system it does not 
have orbitals with vl close to 1/2. To get convergence, either with the DDI 
or the DD2 prescription, the difference of the total energy calculated with 
different L values, must be zero. This is the case in the DDI approach where, 
even with even-L values, the HF energy differences are almost exactly the 
same as the pairing ones. As a result, both cancel out and the total energy 
is convergent, as expected. With the DD2 prescription and even-L values. 
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Table 2 

Differences AE (in keV) between the projected energies calculated with L = 8 
(Ef) and L = 20 (Ef^) in the Fomcnko expansion for the nucleus -'^^''Er taking into 
account all the contributions of the force to the fields. [DDI] and [DD2] are meaning 
the prescription used for the density-dependent term. In the last column ^Ef^ is 
the difference between the projected energies calculated with L = 9 and L = 21. 

however, the HF results are different from the pairing ones and as a result no 
convergence for the total energy is found. We know, however, that the DD2 
divergence is integrable, that means if the integral is carefully performed the 
convergence must be found. This is what happens if we do the calculations 
with odd-L values (see last column) where the largest deviation found is 3 
keV. Obviously, more sophisticated integration methods will provide better 
convergence. 

An interesting point concerning the density prescriptions DDI and DD2 is to 
know if there is a big difference in the results calculated with DDI or DD2. 
In Fig. 4 we show the 7-ray energies along the yrast line of the nucleus ^^^Er 
in the HFBc with the DDI (in the figure PAVl) and with DD2 (in the figure 
PAV2). As we can see the behavior is very similar and the physics is clearly the 
same. A similar conclusion was found in [15] in the Lipkin-Nogami approach. 
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Fig. 4. E^[MeV] versus angular momentum for both prescriptions of the density. 

3 The Vciriation after projection method : The nuclei ^^Cr and 
50Cr 



In ref. [24] and [25] a good agreement between the HFB approximation, the 
exact shell model diagonalization and the experimental results was found for 
the rotational states of the nuclei '^^Cr and ^°Cr. Only the 7-ray energies 
along the yrast bands in the HFB approach were systematically around 0.5 
MeV smaller than experimentally observed. A simulation with the exact shell 
model diagonalization in [24] showed that the shift was due to a deficient 
treatment of the pairing correlations of the HFB method in the weak pairing 
regime. It is well known that the variation after projection method is able to 
provide pairing correlations even at the limiting case when the HFB approach 
collapses to the HF one. In this section we shall investigate the mentioned 
nuclei in the VAP method. 
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Fig. 5. Pairing correlation energies in the HFBs and VAPg approaches for the nucleus 
^^Cr and ^°Cr respectively. 



The HFB calculations of refs. [24] and [25] were performed in the standard 
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approximation, i.e. HFBs. For nuclei like these with little pairing -low level 
density at the Fermi surface- is very unlikely to find orbitals with vf, = u\, 
nevertheless we have checked during the minimization process of the VAP^ 
that no orbitals were populated in that way. Additionally, we have performed 
also the HFBe and the VAPg calculations. In Fig. 5, we represent the HFB^ 
and VAPs (see eq.(B.14)) pairing energies for the nuclei ^^Cr, left hand side, 
and ^'^Cr, right hand side. In the HFBs approach the pairing energies are rather 
small, for *^^Cr practically zero for all spin values, with the exception of the 
7 = and 2h values where they are very small. For ^°Cr they are also very 
small and only for 7 = 0, 2, 4 and VIU differ from zero. In both nuclei we find 
the Coriolis antipairing effect that leads to a sharp transition from a (weakly) 
superfiuid system to a normal one. The VAPs solutions, on the other hand, 
are well paired and no sharp transition is seen, only a smooth decrease in the 
correlation energy is found. 




1.0 2.0 3.0 1.0 2.0 3.0 4.0 5.0 
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Fig. 6. The transition 7-ray energies versus the angular momentum in the HFE^ 
and VAPs approaches for the nucleus ^^Cr and ^°Cr. The experimental results are 
also shown. 

In Fig. 6, the 7-ray energies along the yrast band are displayed in the HFBs 
and VAP5 approaches together with the experimental results for ^^Cr, [26] 
and ^°Cr, [27]. As discussed in [24], the HFBs 7-ray energies for ^^Cr show 
the same trend as the experimental ones but are about 0.5 MeV smaller. The 
VAP5 approach, on the other hand, due to the larger pairing energy leads to 
smaller moment of inertia improving the agreement with the experiment. The 
VAPs results come closer to the experiment and the backbending is slightly 
better described. For "'''^Cr the situation is similar, the VAP^ approach improves 
considerably the HFB^ one specially at low and medium spins. One may ask, 
however, why we do not get enough pairing correlations in the VAP^ method 
as to match the experimental results. There are several possible answers, first 
in our calculations we do not include p-n pairing which is known [28] to play 
an important role in these nuclei. Second, the fact that in these nuclei, in the 
HFB approximation, there are barely any pairing correlations indicates that 
the potential energy curve is rather flat against pairing fluctuations. This kind 
of fluctuations is not included in the PNP method. 
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In Figs. 7 and 8, we finally present the results of the exact HFB calculations, 
HFBe, and the corresponding VAPe method. The pairing correlation energies 
are represented in Fig. 7. As expected [18], the neglected pairing terms in 
the HFB^ approach have an antipairing efi^ect, mainly the Coulomb term. 
The small pairing energy of the HFE^ approach is quenched in the HFEg, in 
particular the proton pairing energies completely vanish. The pairing energies 
in the VAPg approach display this effect even more clearly. The total VAPg 
pairing energy is considerably reduced with respect to the VAP^ approach. 
Furthermore, while in the VAPg the proton pairing energy is always larger 
than the neutron one in the VAPg the proton pairing energy becomes always 
smaller than the neutron one. 
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Fig. 7. Pairing correlation energies in the HFBe and VAPg approaches for the nuclei 
^®Cr and ^'^Cr respectively. 

In Fig. 8 we finally present the 7-ray energies in the HFBg and the VAPg 
together with the experimental ones. Since the pairing energies are negligible 
in the HFBg and HFB^ for both nuclei, both sets of 7-ray energies practically 
coincide. With respect to the VAPg values they are of the same quahty as the 
VAPs though somewhat smaller than these one, as expected from the behavior 
of the pairing correlations. 



7 Conclusions. 



In this paper we have derived the expressions needed to perform a particle 
number projected HFB calculation with a finite-range and density dependent 
force in the triaxial basis. We have studied and thoroughly discussed the di- 
vergences that appear for special values of orbital occupancies when exchange 
terms are neglected. We have performed exact and approximate calculations 
(neglecting some exchange terms) and found that in the latter case some prob- 
lems associated with the convergence of the particle number projection arise 
leading in some cases to an inaccuracy of about 10 % in the values of the 
7-ray energies. Interestingly, all terms of the force, even of the Coulomb part. 
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Fig. 8. The transition 7-ray energies versus the angular momentum in the HFBg 
and VAPe approaches for the nuclei '^'^Cr and ^'^Cr. The experimental results are 
also shown. 

must be taken into account. This feature presents a challenge for interactions 
where a given force is used in the particle-hole channel and a different one in 
the particle-particle channel. 

We have furthermore performed variation after the projection for the ^'^'^Od- 
nuclei with the Gogny force. We found strong pairing correlations, at vari- 
ance with the weekly correlated HFB approach. These correlations smoothly 
decrease with increasing angular momentum indicating that no sharp phase 
transition takes place as predicted by the HFB approach. The projected calcu- 
lations improve the HFB ones but some correlations (probably T=0 pairing) 
not present in our approach are still missing. 
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A Appendix A: The Gogny Force. 

In the calculations we use the Gogny interaction [29] as the effective force . 
The main ingredients of this force are the phenomenological density dependent 
term which was introduced to simulate the effect of a G-matrix interaction and 
the finite range of the force which allows to obtain the Pairing and Hartrce- 
Fock fields from the same interaction. We use the parametrization DIS, which 
was fixed by Berger et al. [30]. The force is given by 
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vi2 = E e-(^^-^^^)'/'^'(H^. + B,P^ - H,Pr - M,P„Pr) + 

i=l 

+ Wls{^i + ^2)k X (5(fi - f2)A; + 

+ h{l+x^P^)5{n-r2)p^'^ (^(n +f2)) , (A.l) 
and the Coulomb force 

vg = {l + 2n,){l + 2r2,)- — . (A.2) 

The density operator, /3(f^, is given by 

p(f) = x:5(f - ri) = 5:0*(f)0,(f)(5.i5,)cic, = 5:/,,.(f)cic, . (a.s) 

1=1 ij 

In the two-body interaction used in the calculations we also include the one- 
body and two-body center of mass corrections. 




(A.4) 



B Appendix B: Calculation of the projected energy in the canon- 
ical basis. 

In this appendix we show that the projected energy is always well-defined if 
the direct (Hartree), exchange (Fock) and pairing terms of each component of 
the force are taken into account, i.e. if no exchange terms are neglected. 

The projected energy can be written in the following way Q 

E^ = Tr{tp^) + E^,f + E^ (B.l) 

where 

k\k2kzki I 

K=\ E ^krk2ksk,Y.yi'^klk2i'^i)'^kik,i'^i) (B-3) 

kik2k'^k4 I 

^ To simplify the expression we consider the projection onto good particle number 
only in one dimension in the gauge space. 
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These energies are given in terms of the (/j^-dependent density matrix p{(pi) 
and the pairing tensors K^^{ipi) and K^^{ipi). These quantities can be easily 
evaluated in the canonical basis. In this basis, denoted by {a, a''^}, the usual 
density matrix pfcjfcj = {^\a\^ak^\^) is diagonal and the usual pairing tensor 
i^kik2 — (^lofcaOfeil^) has the canonical form. In this basis the quantities we 
are interested in are given by 



Pkik2W) - 



($1 


l4i 









These matrix elements can be easily evaluated 



k>0 

= Y[{-\{uk + Vkakak){uk + Vke^"^'alal)\-) (B.b] 

k>0 



m4+ 



2 2i^l N 



fc>0 



here k denotes the canonical conjugated state of k. To calculate {^lal^akj^e^'^^'^ \ ^) , 
we shall proceed in the same way. Let us first assume k2 ^ ki, then 

($|4^afc,e^'^'^|<l>) = l[H{uk + Vka-kak)alak,{uk + Vke^''^^alal)\-) = 0,(B.6) 

fe>0 



that means, the (/9;-dependent density matrix pk-^k,^{Lpi) , like the normal density 
matrix pfc^fcj, is diagonal in the canonical basis. For k2 = ki = k we obtain 



($|4a,e'v''^|$) = ^^ e^^^^' ■ J] + ^^e^^'^O- (B.7) 

m>0,m^fc 

Therefore, in the canonical basis, the matrix pkm{vi) is given by 



^2 . ^2m 

Pkrai^l) = hm " ^2 ^ ^2^2^ ' (^-8) 

We see in this expression that for Uk — Vk and (pi — 7r/2, the matrix element 
diverges and one has to be very careful in the calculations. It is easy to show 

that pkrniVl) = Pkmi^l) ■ 

Let us now calculate the (/^rdependent pairing tensors in the canonical basis 
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in the same way 



($|afc,a,,e'^'^|$) = S^^^^u^^vu, IJ (^^ + ' e'^^O- (B-IO) 

m>0,m7^fe2 



The non-zero matrix elements are given by 



obviously, = "'^ijfc ^^<i '^fel — "'^^r Again the dangerous terms are for 
Uk = Vk and = 7r/2. 

^From eq. (B.2) and using the expressions derived above, we obtain 



^hf — 2 5Z '"kik^kski X] yiPkki'Plj^kskhikPk'k'i'Pljhik'hik' 

kik2k3k4 I 

= 2 S 'YlyiPkk{fl)Pk'k'{fl) 

= ^ (^^fefe'fefc' + %fe'feik' + %fe'fefe' + Vkk'kk')^yiPkk{vi)pk'k'{(Pij'B.12) 

^ k,k'>0 I 

Taking into account that yi{<fi) ~ ($|e"^'^|$), we see that the possible poles of 
Pkki'^i) and Pk'k'i^i) are canceled when k ^ k', i.e. we may only have problems 
if A; = k'. The contribution to the Hartree- Fock energy in this case is given 
by 



[Ehf]pole = E E ^kkkk yi , 2 , „2^2i^ 

fe>o I y^k + 



=A^EE^.^.^7rra^- n iul + vl^e^'-^). (B.13) 

fe>0 Z y^k + ^fe«^ "^'yl m>0,m7^fe 



^1 

. This con- 



where we have made use ofeq. (5) and AT = Ef=i('^'|e'^''^"» 
tribution clearly diverges for ul = and ipi = n/2. The pairing energy term 
is given by eq. (B.3), in the canonical basis it takes the form 
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^pai ~ 7 X! X! yii'^kkk'k''^'ik'^k'k' + '"kkk'k''^k\'^h'k' + 
k,k'>0 I 

+ ^kkk'k''^kk'^k^k' + ''^kkk'k'f^ki'^k^k') 

= E Y.yi'^kkk'k'i^WA'- (B-14) 

fcfc'>0 i 

As before, cancels one possible pole, and ii k — k' the contribution to the 
energy is 



V^g \pole = Z^Z^ ^kkkk Vly^ o o\,r. \2 

=ArEE--.u.7-^tk^- n + (B.15) 

fc>0 / \^k^ ^k*^ ) m>0,m^k 

which also diverges for u\ = vl and ipi = ti/2. The total contribution to the 
energy is given by the sum of eq.(B.13) and eq.(B.15), we obtain 



[^hf]pole + [Eg]pole — '^kkkk ^ ^ 



k>o I yul + vle^-^'P') 

MT.H^kkk-kvle'''''- n i^n + vl^"'''') (B-16) 

fe>0 I m>Q,m^k 



which clearly does not diverge at ul = vl and = 7r/2. 

This result shows that in case of particle number projection one cannot arbi- 
trarily neglect exchange terms of any component of the two-body potential, 
including the Coulomb one. 



B.l Density-dependent interactions. 



In the previous demonstration we have not considered the case of density- 
dependent interactions v^^k^k^k^ in the analysis of the convergence. Wc shall 
distinguish between the two choices for the Hamiltonian density in order to 
study this problem. 



B.1.1 Projected density. 

Taking into account the expression (30), one sees that the density-dependent 
energy depends only on the proton and neutron projected densities. We must 
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therefore analyze the behavior of these densities in the presence of poles. The 
projected density is given by 

Pk,k2^12yiPkik2{Vi)- (B.17) 

As we saw in the previous section, any pole of pk^k'^i'^i) will be canceled by a 
zero of the same order from the factor yi. Therefore, if we take the projected 
density for the density dependent term, the total energy including this term 
is well-defined. 

B.1.2 Mixed density. 

In this case the energy of the density-dependent term is given by eq. (34). 
As we can see in this expression and in eq. (36), the densities p{(pi^) and the 
norms yi^ appear pairwise indicating a cancellation of possible poles and zeros. 
However, problems may arise if the interaction Vdd itself has some poles. In 
this prescription Vdd depends on p(f, ipi) given by eq. (33). Writing the density 
matrix elements in the canonical basis we find that 

P{r, ^i) = E ih^ (rl + f-k-k (r-)) (B.18) 

k % + ^k^ 

where the functions fkk'if) are defined in eq. (A. 3). Clearly, if one of the v^, 
let's say v^^, equals 1/2 then the density is singular at ipi = tt/2 and in the 
neighborhood of this point it behaves as 

P(r, Vi) - -\ {fk^ko (r-) + h-,^ (r)) (^^^^ (B.19) 

which is clearly divergent. This singularity, however, does not pose any real 
drawback as the density dependence of the interaction is proportional to 
p(f, ipi)°' with a = 1/3 and therefore the singularity is integrable. 
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